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Abstract. The paper proposes a method to obtain the optimal basis set for solving 
the self consistent field (SCF) equations for large atomic systems in order to calculate 
the energy barriers in tunneling structures, with higher accuracy and speed. Taking 
into account the stochastic-like nature of the samples of all the involved wave functions 
for many body problems, a statistical optimization is made by considering the 
covariance matrix of these samples. An eigenvalues system is obtained and solved 
for the optimal basis set and by inspecting the rapidly decreasing eigenvalues one may 
seriously reduce the necessary number of vectors that insures an imposed precision. 
This leads to a potentially significant improvement in the speed of the SCF calculations 
and accuracy, as the statistical properties of a large number of wave functions in an 
large spatial domain may be considered. The eigenvalue problem has to be solved only 
few times, so that the amount of time added may be much smaller that the overall 
iterating SCF calculations. 

A simple implementation of the method is presented for a situation where the 
analytical solution is known, and the results are encouraging. 
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1. Introduction 

Ab initio methods become more and more efficient in various scientific and technical 
apphcations, as the involved physical principles, numerical methods and computer 
hardware display a constantly improvement. The whole effort is sustained by the big 
promises of these methods in the general field of the computer simulation of matter 
properties, with applications in physics, chemistry, biology, and many interdisciplinary 
scientific researches. However, there are still many problems that demand better 
solutions at any level of the ab initio methods and their solving is limited by the 
enormous computational effort implied by these kind of calculations even for the 
very small clusters of atoms that can be dealt today. As the large accessibility of 
supercomputers will be probably delayed for an unknown period and since there is a 
continuous need for larger atomic systems calculations, some progress is certainly needed 
in both physical principles (more accurate and elaborate models) and numerical methods 
(more precise and fast algorithms) to achieve this goal. 

The physical part of this scenario has already an eight decade history, starting with 
the birth of the quantum mechanics. The Hartree method of the self consistent field 
(SCF) founded in the third decade of the last century was soon amended to include 
the exchange integral leading to the Hartree-Fock (HP) equations which still remain 
the least empirical way for ab initio calculations. Their well known expressions reveal 
the necessity of an iterating process, as each wave function depend on the others and 
they are present both under the differentiation and integration operators. Considering 
the Born Oppenheimer approximation, the left side of the equation for each particle is 
composed by the one electron term, the coulombian interaction between electrons and 
the exchange term due to spin: 



where (r) is the one electron wave function, U (r) is the potential energy in the nucleus 
field, r and r' are the vectors to the nucleus of the two electrons interacting and Sj is 
the spin index. As usually, for the sake of simplicity, natural system units have been 
used in this equation 

About three decades have been necessary for providing a solid theoretical base, 
in the framework of the Kohn-Sham (KS) formalism [Ij of the Density Functional 
Theory (DFT) , for including the correlation term, and other decades for its satisfactory 
evaluation. Using the properties of an homogenous electron gas and introducing the 

N 

functional of density, for a sufficiently slow varying density n (r) = ^ |\E'i(r)| , the 
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system of equations reads [T] : 

j-iv^ + ^ (r) + /i, (r)| VP, (r) - J dr'"^^^^^, (r') = e,^, (r) (1.2) 

where ip (r) includes the one and two electrons potential, /ic (r) includes the correlation 
effects and the last term of the left side of the equation is the new form of the exchange 
correction. Although theoretically more accurate than HF, the KS method has some 
necessary simplifying assumptions of the new parts of the model so that it is often 
considered slightly empirical in these aspects. For the most situations this formalism 
is highly efficient and the theoretical improved accuracy is present in many of the 
applications. However, probably due to its empirical part, it still fails sometimes, as HF 
also does in other situations due to the neglecting of the correlations. 

Finally, it is worth to mention the existence of many other more or less empirical 
methods, with a greater speed, but with a more limited scope, with codes commercially 
available and widely used in various applications. As the physical model is still a rather 
simple one, refining the model is not expected to improve the speed of the ab initio 
calculations, but rather their accuracy. 

Many other post HF models have appeared in the last decades, offering a multitude 
of choices for various aspects of the ab initio calculus. Thus, the exchange term may have 
several forms, as: exact HF exchange. Slater local exchange functional [2], Becke's 1988 
non-local gradient correction to exchange [3], Perdew-Wang 1991 generalized gradient 
approximation non-local exchange [1],[5]. For the correlation term the most accurate 
expressions seems to be: Vosko-Wilk-Nusair (VWN) local correlation functional [6], 
Perdew and Zunger's 1981 local correlation functional [7], Lee- Yang-Parr non-local 
correlation functional |8j|, Perdew-Wang 1991 local correlation functional |4j. 

Concerning the numerical methods, the algorithms and the mathematics used in 
the ab initio calculations, there are many popular techniques that may be considered, 
each of them with their pros and cons [9] . Here the possibilities are more diversified, as 
the form of the self consistent equations may be purely differential or integro-differential, 
due to the exchange and correlation terms. 

The differential form is often treated by the Numerov's fifth order method, which is 
robust and accurate but is not self starting and require some initial iterations, as many 
other point by point methods of high order. A notable exception should be the forth 
order Runge-Kutta method but it is not well suited for boundary conditions equations 
as the HF ones. Some shooting method must accompany the point by point methods 
and, although this provides the eigenvalue of the equation (which sometimes is the main 
goal), it imphes an iterating process that leads to a huge amount of computing effort. 
Furthermore, it must be used both for the wave function equation and for Poisson 
equation for finding the Hartree-Fock potential generated by the charge density. 

An important step for improving the ab initio methods' efficiency was made by 
Roothaan [10] who transferred the calculations to linear algebra in the form of a 
generalized eigenvalue problem, using non orthogonal basis set. The HF equations 
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reduce to the following matrix equation which is more suitable for a numerical 
calculations: 



where F is the Fock matrix, C is a matrix of coefficients of the two electrons interactions, 
S is the overlap matrix of the basis functions, and e is the matrix of orbital energies. 
Thus, a class of new numerical techniques became eligible and an increasing interest for 
appropriate choosing of the basis set emerged. 

Among the various methods with a reduced need of iterations that are currently 
used there are: finite difference method, finite element method, Galerkin method, 
collocation method, etc. The most promising class of methods for such boundary 
conditions equations is considered to be the class of spectral and pseudospectral methods 
[12], as they have already been successfully used in various other fields. Their 
evanescence property (exponential decay of the error with the number of sampling 
points) and the absence of iteration processes are very attractive, but the main problem 
is the bad conditioning that often appears in the linear algebra implied. However, the 
matrix conditioning numbers are generally satisfactory if a low dimension of the vectorial 
space is chosen, but this tends to increase the errors if the basis set is not an optimal 
one. 

We consider that the spectral methods in the ab initio calculations have not been 
entirely exploited, as even the basis sets seem to be still "empirically" chosen: by 
subjective considerations (as Chebyshev polynomials are the usual basis of choice |llj). 
or artificial approximations that facilitate the analytical calculations with a great impact 
on the speed of the numerical part (for example the Gaussian basis) . As the number of 
vectors used by a basis must be always finite and as small as possible to minimize the 
matrix condition number and the overall computational effort, the problem seems to be 
the choice of an optimum basis that provide the minimal error when the dimension of 
the vectorial space of the wave function representation is highly reduced. 

In the next chapter we propose a method for properly choosing the basis set for 
achieving this minimization of the errors when dealing with a smaller number of vectors 
that is usually needed. An example for the implementation of the method will be later 
presented and some conclusions and further suggestions will be drawn. 

2. The stochastic nature of the wave functions' samples for many body 



The radial wave function in a hydrogenoid atom, which is the primary natural 
approximation of the many body wave function has the well known form |13j : 



where n is the principal quantum number, £ is the angular momentum number, r is the 
distance to the nucleus in Bohr radius and (r) are the associate Laguerre polynomials. 



FC = SCe 



problems 




(2.1) 
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Figure 1. The various wave functions in a hydrogen- like atom for n = 1, 2, 7 

Every electron moving in the potential determined by the charge distribution 
created by such type of wave functions is thus influenced by all the others electrons. 
For systems with many electrons, the family of functions has the appearance shown 
in figure [1] (generated using Mathematica software). The asymptotic behavior of the 
wave functions allows one to deal with a finite interval, but this interval seems to be 
quite large for the usual methods for point to point differential equations . Thus, the 
errors are very important for most of the orbitals that imply distances above 15-20 Bohr 
radius, but the influence of the further regions is clearly present in the real system, and 
needs to be considered. 

If one samples the distance to the nucleus, a pseudo-random distribution of the 
local values of these wave function appears as shown in figure El 

Although the entropy displayed by these values is not very high, it is possible to 
consider a stochastic influence between the electrons, and use the second order statistic's 
methods for dealing with them. The concept of the "mean field", with a clear stochastic 
connotation is present in many of the physical models accepted now-a-days but this 
context has not been considered enough in the present theories, except for some low 
order Monte-Carlo methods. 

3. Using the Karhunen-Loeve theorem for calculating the optimum basis 
set 

The spectral methods used for solving the eigenvalue problems in simple or generalized 
form, first approximate the solution as a linear combination of continuous functions - 
basis vectors- and then plug in this solution in the original equations to determine the 
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Figure 2. The equidistant samples of the wave functions in a hydrogen-like atom for 
n = 1,2,...,7 



coefficient matrix. 

N 

y{x) ^^dcpiix) (3.1) 

The domain of the independent variable is sampled in a number of points xj (as 
in the collocation method) where the equations are supposed to be satisfied, and the 
resulting linear system is solved for the coefficients. 

N 

yi^j) = '^Ciipi{xj), j = 0,...,N (3.2) 

1=0 

As the values of y (x) are known for the boundary points, a system of + 1 
equations is obtained, and its solution gives the unknown coefficients. Although it seems 
to be a very simple procedure, there are several problems that may occur and must be 
somehow taken into account. First of all, if one uses a simple, equidistant sampling, 
the convergence of the method, theoretically exponential, is lost due to the Runge 
phenomenon. That is why various non equidistant sampling methods are currently used 
with a smaller step at the extremities of the interval (as Chebyshev, Gauss-Lobatto 
or Legendre points). The second, and the most serious problem is the fact that the 
relation ( 13. ip is still an approximation and eventually the equality is true for a huge 
if not an infinite number of basis vectors, for the most part of the applications. For 
practical purposes, only a very limited number of these basis vectors may be used, for 
two reasons: 

- The computational effort increases with at least A^^, limiting the number of wave 
functions that may be dealt with reasonable timing. 

-Increasing A^ always increases the condition number of the matrix and hence serious 
errors occur for A^ values above 10 - 20. 
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That is why the basis set must be very carefully chosen, as it may prove itself as 
crucial for the performance of the whole calculus. 

An important observation must be now taken into account: the coefficients of the 
linear combination (13.11) carry some amount of redundancy in all but one case. This 
unique case implies that all of them to be uncorrelated, in the second order statistics 
meaning: 

E {dCj) = \iSij, i,j = 0, (3.3) 

where the symbol E (.) stands for the expected value operator and Sij is the Kronecker 
symbol. 

It follows that if the coefficients are uncorrelated, the information lost by truncating 
the series (13. ip to a number M < N oi terms (for the above presented reasons) is 
minimized, and the truncating errors are also minimized. 

Hence, a basis set that ensures a pairwise uncorrelated set of coefficients must 
be found, which is possible by using the Karhunen-Loeve theorem [T3] [15]. It states 
that the basis set that ensures the minimization of the errors due to the truncation of 
the decomposition expression of a continuous function is the solution of the integral 
equation: 

j Kyy (x, x') Lfi (x') dx' = \iipi {x) (3.4) 

D 

which is a second kind homogeneous Fredholm equation, i.e. an eigenvalue problem. 
Here D is the orthogonality domain of the eigenvectors ^pi {x) that form the optimal 
basis set, Aj are the associate (positive) eigenvalues and Ky (x, x') is the autocorrelation 
function of the initial function y (x) defined by: 

Kyy{x,x') = E{y{x)y{x')) (3.5) 

The set of orthogonal functions obtained from (13.41) is complete and the functions are 
square- integrable. 

If the initial function y [x) is the wave function , one may consider that the 
eigenvalues represent the probability density associated with each mode, and according 
to the KL theorem the set of basis functions is the optimal from all possible sets, i.e. 
the decomposition series converges as rapidly as possible. The method also minimizes 
the representation entropy and is equivalent with the minimization of the mean-square 
error resulted from the truncation. 

It is widely applied in signal theory (detection, estimation, pattern recognition, 
noise rejection, data compression for storage and image processing), physics (stochastic 
turbulence processes [H]) and biology [17], under various names: Uncorrelated 
Coefficients Series, Principal Component Analysis |18j. Hotelling Analysis |18j . 
Quasiharmonic Modes [17], Proper Orthogonal Decomposition [16], Singular Value 
Decomposition [19] etc. 

The stochastic-like nature of the samples of all the wave function in a many body 
problem, suggests that this method could be also applied to the SCF problems. Since 
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the numerical processing involves the discretization of the implied functions, the KL 
transform is often met in matrix form with finite dimensional vectors. Considering a 
total number of N^j hydrogen-like wave functions and sampling all of them in A^^ points, 
we obtain column matrices each of them containing Ng elements: 



J = 1,2, 



AT,, 



(3.6) 



The resulting A^^ x matrix 
Y=(Yo Yi . . 



(3.7) 



contains all the samples of all the wave function for hydrogen-like atoms. By averaging 
the similar samples and subtracting the result from each sample of each wave function 
one obtains a matrix with columns of zero centered values, as the KL transform demands. 



(YJ 



2 



Yat^) 



with 



(Y^ 



Vj [Xi) 



Xi 



0,1,..., AT,; j 



0,l,...,Ar„ 



(3. 



(3.9) 



Using the centered samples matrix Y*^ and its transpose [Y*^]"^ one may construct 



the Ng X Ng covariance matrix defined by 



K 



yy 



E 



(3.10) 



Using these matrices and sampling the independent variables x and x' in eq. (13. 4p 
it is transformed in an eigenvalue system for the covariance matrix: 



K $ 



1,...,N, 



(3.11) 



which has to be solved to obtain the column matrices $j containing the samples of 
the needed basis function ipi{x) , i = 0,1,. ..,N . One should take sufficient number 
samples to assure A^^ > A^ + 1 , where A^ is the dimension of the primary space 
of the decomposition. The eigenvalues Xj obtained by solving the covariance matrix 
are supposed to be strongly decreasing with j , and this was numerically checked, as 
shown in the next chapter. The first M largest eigenvalues indicate the corresponding 
eigenvectors ipi , with i = 0, 1, ...,M that should be taken into account in a truncated 
decomposition with acceptable errors. The more rapidly the eigenvalues decrease, the 
smaller is the number of vectors in the basis set that must be taken into account. In the 
next chapter we will show a simple example of implementation of this technique, and 
some preliminary results for the covariance matrix and the basis set. 
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Mathematically, the KL equation fl3.4l) is equivalent to a transformation which 
diagonalizes a given matrix K and turns it to a canonical form 

Kyy = VAV (3.12) 

where A is a diagonal matrix. Indeed, by forming a matrix $ , who's columns are the 
basis vectors $j 

$ = ($1 $2 . . . $Nw) (3.13) 

we may define a transformation from the primary space of the samples yj (xj) to a 
secondary space as: 

Z = (3.14) 

In this space the covariance matrix defined similarly with eq. (I3.10p has a diagonal 
form: 

Kzz = E (ZZ^) = E ($^Y {^^Yf) = ^^E (YY^) * 



(3.15) 

= #^Kyy$ =Diag(Ai) 
where we took into account eqs. (13.111) . fl3.14p and the orthogonality of the basis set 
= I. It then makes possible to use one of the many well developed methods for 
diagonalizing the covariance matrix in order to find the eigenvectors and to select those 
of them which correspond to the largest eigenvalues as the optimal basis set. 

4. Example of implementation of the method and some results 

In order to test the above described procedure we considered a simple situation for the 
ground state hydrogen atom which has an analytical solution for the wave function and 
a precisely known orbital energy. This allows an objective test by comparing our results 
with an exact one and is also a rather difficult situation, because the covariance matrix 
is constructed using all the 28 different orbitals for = 1 — 7 and i = — 6. Thus, the 
IS orbital is at the extremity of the whole range of wave functions and better results 
may be expected for an intermediate orbital. The radial part of the wave function may 
be obtained analytically by imposing the boundary conditions y{x) = i/i = for x = 
and y{x) = y/ for x = b {in Bohr radius) to the radial part of the Schrodinger equation 



^y" (^) + 



Z ^ £(£+!] 



^ '0 2^ yix) = Eny{x) (4.1) 
X + e 2x^ + e 

where x = r/ag, y{r) = rR{r), oq is the Bohr radius and, in order to avoid the singularity 
in the origin in the numerical calculations, we added the very small constant e (say 10^^'^) 
at the denominator of the coulombian and orbital terms (that does not affect seriously 
the numerical results). 

For K-shell electrons and Z = 1, the solution in a the region (0, b) with b finite, may 
be expressed in terms of exponential integral functions Ei (z) as the symbolic software 
generated expression: 

_ e^-^ {e2(^+^)e -(x + e) [e^' - 2eEi (2e) + 2£Ei (2x + 2^)] } 
^ ^""^ " e2(6+^)£ + 2eib + e) Ei (2^) - (6 + e) [e^' + 2eEi {2b + 2^)] ^^'^^ 
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Figure 3. An example of the obtained covariancc matrix Kyy for A^^ = 6 samples 
and = 28 wave functions 

0.999925 -0.0121609 -0.0003S27S6 -0.00108047 -0.000781293 -0.00030386 

-0.0119481 -0.983512 -0.171783 0.0469124 0.0291472 0.00113039 

0.000504555 0.118928 -0.841887 -0.509929 -0.128696 0.0221725 

-0.0024028 -0.111528 0.385493 -0.516731 -0.645344 -0.394299 

-0.00103903 -0.0762847 0.324052 -0.58617 0.251859 0.694364 

0.000133316 0.0124092 -0.0900129 0.3566 -0.709003 0.601576 



Figure 4. An example of the KL transform matrix for Ng — 6 samples and 
= 28 wave functions. Each column contains the samples of a function $i of the 
optimal basis set 

{4.31474, 0.083716, 0.0186048, 0.00659303, 0.00281989, 0.000764637, 0.000170135, 9.34579x10-', 
3.56925x10-'', 1.68596 x10-% 1.62426 x10-^", 4. 91027 x 10-^% 4.92398x 10"^', 7.54075x 10"^', 
-2.06174x10-^', -3.88258x10-^*, 2.14183x10-^*, 1.19703x10-^*, -9. 75731 x 10-^% -2.83214x 10-^'} 

Figure 5. An example of the eigenvalues Xi obtained for Ns — 20 samples and 
= 28 wave functions 



We used the natural units and thus the energy is expressed in Hartree {E^ = —0.5 
Hartree for Z = l,n = l,i = 0). The boundary condition yf, appearing as a global 
coefficient of the solution may be taken arbitrary or may be determined from the 
normalization of the resulted wave function. 

Using a parameter 6 = 40 Bohr radius and the samples of the 28 wave function we 
construct the covariance matrix according to eq. (13.101) . We exemplify this matrix only 
for 6 equidistant samples for each wave function, for space saving in figure [3l 

By solving numerically the eigenvalue problem for this covariance matrix we obtain 
the KL transformation matrix $ presented in figure HJ also only for 6 samples per wave 
function. 

For a realistic calculation, we used Ng = 20 equidistant samples and from the 
covariance matrix (not shown here) we obtained the eigenvalues listed in figure [5l 
The equidistant sampling method was chosen for simplicity of this exemplification 
but increasing errors are expected at the extremities of the interval due to the Runge 
phenomenon. Of course, in practice we recommend more appropriate sampling methods 
like those mentioned in chapter 3. 

One may see that the eigenvalues are indeed rapidly decreasing, confirming the 
possibility of reducing the dimension of the basis set, as it was anticipated in the previous 
subsection. It follows from the presented values that only 8-10 eigenvectors should be 
considered, the others corresponding to much smaller eigenvalues. 
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For solving the differential equations using spectral methods, one needs the 
derivatives of the basis function. Taking into account that these function are known 
in a small number of points (their samples) it is not effective to use a direct numerical 
differentiation. Instead, one may interpolate those samples for each of the basis functions 
and thus obtain an analytical expression, for example a polynomial. For the sake 
of simplicity, we used a Lagrange interpolation (but some other methods may be 
investigated for better results) and obtained the basis functions as exemplified in figure 

El 

Using these analytical form for the basis set, we checked the numerical solving of 
the radial differential equation 14.11 with the following parameters: i = 0, n = 1, En = 
—0.5, A^s = 8, a = 0, b = 7, yf = 10"'', e = 10^^°. We obtained a wave function as 
presented in figure [3, which is very close to the analytical one presented dashed in the 
same figure. 

The errors introduced by the interpolated basis set, obtained as the difference of 
the two members of the differential equation are presented in figure [HI where one may 
see the effect of the uniform sampling that we used for this example: the errors are 
highly increasing towards the boundaries. Also, the result of the bad conditioning of 
the system (13.21) may be noticed for high values of x both in figure [7] and figure M It is 
possible that this conditioning should be improved if proper sampling and interpolation 
techniques are used. Anyway, the result is very encouraging, taking into account that 
for other basis sets (for example gaussian), in similar conditions, clearly higher errors 
are displayed. 
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Figure 7. A comparison of the analytical solution (dashed) with the numerical one 
(continous) obtained for iV^ = 20 samples and Nyj = 28 wave functions 

e 




Figure 8. Difference between the two members of the radial equation obtained for 
Ns — 20 samples and = 28 wave functions 

5. Conclusions 

Although the above presented method is highly susceptible for improvements, the simple 
example presented for implementing the concept of the optimal basis set obtained using 
the covariance matrix over the set of typical wave functions gives satisfactory results. 
The most important feature of the method is the virtual increase in speed due to the 
theoretically smallest number of basis functions needed and their construction using the 
information contained in a big number of possible orbitals. 

The price paid is the necessity to calculate and diagonalize the covariance matrix, 
but it is important to notice that it must be done only at the start of the iteration 
process and maybe in some intermediate steps. The assertion above is experimentally 
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proved by the fact that in data processing an initial set of samples may be used for 
constructing covariance matrices adequate for other different situations. 

Improvements of the presented example may be made in the sampling and 
interpolations processes and the principle may be also applied to other various spectral 
methods used by ab initio calculations. 
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